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Abstract. Derivation of the following formulae for solar position as seen from orbiting 
planet based on a simplified model: sunrise direction formula, solar declination for¬ 
mula, sunrise equation, daylight duration formula, solar altitude formula, solar azimuth 
formula. Use of notion of effective axial tilt, and Rodrigues Rotation Formula, and re¬ 
flections of the orbital quadrants to reduce the general case to the simpler case of the day 
of the winter solstice. Derivation of equations for solar time to clock time conversion, 
and sunrise, sunset, and solar noon times. Implementation of an analemma calculator. 
Comparison of the sunrise direction formula with 304 point dataset of actual sunrise 
data from Earth, obtaining average accuracy of 1.25°, and estimate of Earth’s axial tilt 
of 23.52° — within 0.1° of the currently accepted value 23.44°. 


1 Introduction 

This article first derives a formula for the sunrise direction for a planet orbiting a central sun, in terms 
of the day of the year and the latitude, under a simplified model. Then, building on the method of proof 
used a number of further formulae relating to the solar position [1], [2] are derived. The approach adopted 
reduces the cases of the four orbital quadrants to a single quadrant and then reduces the latter case to a 
single point (the winter solstice). The simplifying assumptions made throughout are : 

1. The planet is a sphere 

2. The orbit of the planet is an ellipse with the sun at one focus^ 

3. The planet radius is negligible compared with its minimal distance from the sun 

4. Whilst orbiting, the planet rotates about a fixed direction axis through its center 

5. The sun can be approximated as a point source of light 

6. A day of the year is approximated as a single stationary point in the orbit at which the planet completes 
a full 360° revolution 


The sunrise direction formula gives the sunrise direction 9 as an angle north of east : 

^ . f sin a cos'0 

9 = — arcsm -;— 

\ COSO 

where : 

9 € [-90°,90°], 

a = planet axial tilt S [0°,180°], 

6 = latitude S (—90°, 90°), 

0 = day of year angle = orbital angle swept out from winter solstice 


( 1 ) 


In a circular orbit/constant speed model the stationary days are evenly spaced around the circle and 
the orbital angle 0 can be defined as ^(360°) for day d G [0, A — 1], where d is the day offset from the 
day of the winter solstice, and N is the number of days in the year. A more accurate relation between the 
day of the year and 0, in the case of Earth, is described in [3] with a formula expressing the angle swept 
out from the spring equinox as a function of time — we can then take this angle at say 12 noon on a day and 
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^The proofs below do not specifically make use of the elliptical property — the only orbital parameter they require is the 
orbital angle ip, thus they would apply to any planar orbit shape so long as the other simplifying assumptions applied. In the 
case of Earth the gravity of the moon and other planets perturbs its orbit away from a perfect ellipse. 
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add 90° to obtain ip, the time in the formula being measured in units of mean solar days of 24 hours duration from UTC 
midnight on 1®* January 2013. In the comparison with real data from Earth in Appendices A.l and B the former model 
is used. To use the latter model instead, the program code which computes the formula’s predictions can be changed as 
described in Appendix C. The former model is a reasonable approximation for Earth whose orbit is very close to circular, 
but for a planet with more eccentric orbit a procedure such as in [3] would be required. 


In the above simplified model the sunset position is symmetrically opposite to the sunrise position, on the west side 
of the horizon, at angle 6 north of west. A negative 6 means south of east/west. An alternative expression is : 


where 


9 = — arcsin 


smai 
cos (5 


sin ai = sin a cos ip 


( 2 ) 

(3) 


defines an ‘effective axial tilt’ ai at orbital angle ip (Figure 5), which makes the geometric situation identical to the 
winter solstice day of '0 = 0°, when viewed from a different angle. The effective axial tilt equals the actual axial tilt a on 
the winter solstice, and falls to zero at the spring equinox, so ai G [0,a]. Any day of the first quadrant QI of the orbit 
is thus equivalent to the winter solstice of Figures I and 2 (or Figures 3 and 4 in the case of the southern hemisphere) 
with the effective axial tilt ai replacing a. Days in (52, QZ, QA are then obtained by simple reflections from Ql. Thus 
with the substitution of ai for a the simpler case of the winter solstice in Figure I applies to the entire orbit, with the 
sun being placed on the opposite side of that in Figure I in the cases of Q2 and Q3. The effective axial tilt is defined 
only for Ql and is in the range [0, a] but if the definition (3) is extended throughout the whole year to a value ai in the 
range {—a, a] then it equals minus the solar declination A (§3). Note the notion of effective axial tilt is only applicable to 
a stationary day at a fixed point in the orbit, since otherwise it would imply the wrong direction of motion for the planet 
in its orbit during the day. 


The formula (I) can also be written as : 

sin 9 cos J + sin a cos ip = D (4) 

or, to estimate the axial tilt from the sunrise direction (Appendix B) : 

sin 6* cos (5 , , 

sina =-;— (5) 

cosip 

provided cosip ^ 0, ie we are not at an equinox, where the formula becomes 0/0. Similarly (4) can be rearranged to 
express day of the year in terms of sunrise position at a known latitude (or Pole Star elevation, Appendix E) if the axial 
tilt is known. 


The term ‘winter solstice’ refers to the standard notion of ‘winter’ of the northern hemisphere, even though towards 
the equator the terms winter and summer lose their normal meanings — for example on the equator maximal solar 
radiation is received at the equinoxes and is minimal at both summer and winter solstices [4]. 

Though astronomical terminology is used throughout, the problems of the solar position under the simplified model 
are purely geometric involving a sphere on an elliptical path intersected by parallel rays from the direction of one of the 
path’s foci. For example the point of sunrise/sunset on a planet is simply when these rays become parallel to the tangent 
plane at the point. The desired sunrise/sunset direction is then the direction of the rays wrt true north [5] in the local 
tangent (or horizon) plane. 

In Appendix A.l, the sunrise direction formula with a circular orbit/constant speed model is compared with actual 
data for the case of the Earth, using 304 data points at 8 latitudes for the year 2018-2019 — the actual sunrise directions 
being taken from www.timeanddate. com [6] and the latitudes from Google Maps. Using an Earth axial tilt of 23.44° [7], 
the formula then matches the actual sunrise directions to within an average error of 1.25°, which is a quite good ap¬ 
proximation for the simplified model. Using the formula’s own predicted axial tilt of 23.52° (Appendix B) the error is 1.26°. 

A chief interest of the sunrise direction formula is how in relating the above quantities measurable by naked eye from 
the Earth using only primitive instruments, the hypothesis of a spherical Earth in a circular orbit and rotating about a 
fixed axis is arrived at using only a geometric argument — since the close correlation with the actual data would be highly 
unlikely to occur unless the above simplified model was a good approximation in practice. Thus it is one method from 
which the ancients could have determined the Earth shape and motions from only naked eye observations and geometry, 
without the benefit of modern astronomy. 

The naked eye measurements corresponding to the above parameters are : 
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(A) n, the day offset from the winter solstice, where the winter solstice is determined by the day on which the Sun is 
at its lowest southerly elevation at solar noon (with the convention of a northerly Sun elevation being > 90°), and 
n is counted from that day starting at zero. For southerly latitudes such lowest southerly elevation would actually 
be in the summer as there the Sun reaches into the northern half of the sky and is low in the northern sky in the 
winter — the day offset would thus be counted from the summer solstice when the Sun is at its highest northerly 
elevation at solar noon, which is at the same time as the northern winter solstice. 

(B) 6, the latitude, is obtained from the fixed angle of elevation of the Pole Star at the given location. Where the Pole 
Star is no longer visible a corresponding fixed point in the southern sky would be used to determine latitude. 

(C) 9, the position of the sunrise, is measured wrt due east, where due east is determined wrt true north, and true north 
is determined from the position of the Pole Star projected onto the horizon orthogonally (the ‘azimuthal’ position). 


The latitude 6 being measured as elevation of the Pole Star is based on the hypothesis of a spherical Earth and very 
distant Pole Star (see Appendix E). Measurement (A) may be recorded by some kind of calendrical system marking 
mid-winter and mid-summer. 

Originally ancient astronomers would only have these observed quantities to work with, gathering data over long 
periods of time, and hypothesizing what may be the larger scale structure that caused them [8], [9], [10], [11]. 

The formula would then relate the above quantities with the right choice of axial tilt, which could be hypothesized 
by some other method or be estimated from actual sunrise data using equation (5) (Appendix B). The constancy of the 
rhs of equation (5) for actual data would give evidence for the hypotheses 1-6 above. 

In short the above mathematical relationship between the observable quantities (A) - (C) gives a strong evidence for 
the hypotheses 1-6 with a circular orbit describing, to a good approximation, the Earth shape and motion through 
space — since it would be very difficult for any other configuration to produce the same correlation. 

The derivation below of the sunrise direction formula envisages a prograde (wrt orbital direction) rotation axis of 
tilt up to 90°. For a retrograde rotation {a > 90°) the axial tilt angle (180 — 9)° could be used, and would be clock¬ 
wise, so the sunrise and sunset positions would simply be swapped, and the formula would remain applicable, since 
sin a = sin(180 — a). The ‘north’ side of the orbital plane is defined as the side from which the orbit appears anti¬ 
clockwise. In using (5) to estimate axial tilt, the result chosen is a < 90° or a > 90° according as the rotation is prograde 
or retrograde, the latter being when the sun rises in the west and sets in the east. From §3 onwards we will only consider 
the case of a G [0°, 90°). 

The sunrise direction formula applies to all latitudes 5 except for the poles, at 90° and —90°, where the sunrise/sunset 
is caused by the orbital motion of the planet rather than the daily rotation — eg on Earth at the North Pole the Sun 
rises once per year on the spring equinox and sets once per year on the autumnal equinox, the opposite being true for 
the South Pole, and this type of sunrise takes much longer than a diurnal sunrise, around 30 hours versus a few minutes, 
for the sun’s disc to fully cross the horizon, and the type of discrete sunrise direction we consider here is not applicable 
[12], [13], [14], [15], [16]. 

For non-pole latitudes beyond one of the arctic circles such non-diurnal sunrises and sunsets also occur at certain 
times of the year — when periods of perpetual day or night are entered. Outside those periods the sunrise and sunset 
are diurnal. These changes occur as the effective axial tilt changes so the actual latitude travels above and below the 
effective arctic circle (the latter becoming a point at the pole on the spring equinox). In the sunrise direction formula 
the periods of perpetual day and night correspond with the arcsin argument going out of range, for there is no sunset 
nor sunrise at these times — and this can be used to determine the range of dates for these periods for a given polar 
latitude, as shown in some examples in Appendix A.2. 
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2 Sunrise Direction Formula 


The formula is first proved for the simplest case of the winter solstice (ie ip = 0), from which the summer solstice case 
then readily follows. The equinoxes are readily checked because formula (1) then gives 0 = 0 at all latitudes which is 
clear from the symmetry — the effective axial tilt is zero, ie. vertical, and each latitude has equal lengths of day and 
night, the sunset and sunrise being along the east/west line everywhere (except the poles which are seeing the non-diurnal 
sunset/sunrise). The general case for the first quadrant Q1 is obtained from the winter solstice case by calculating the 
effective axial tilt using the Rodrigues Rotation Formula (Appendix D and [17]). The cases of Q2, QS, and QA are 
obtained from the Q1 case by considering reflections into these quadrants (Figure 5). 

2.1 Winter Solstice 

We require to show : 

^ , / sina\ 

d = — arcsin -- (6) 

\coso / 


2.1.1 Northern Hemisphere 

The winter solstice case {ip = 0) is shown in Figures 1 and 2, with latitude 6 € [0, 90 — a] in the northern hemisphere. 
The case of i5 = 90 — a means points A and B coincide and formula (6) reduces to — arcsin 1 = —90, ie southerly, as also 
seen from Figure 1. The case of (5 = 0 means B coincides with O in Figure 1 and (6) reduces to the required angle of 
9 = —a. 

The east-west line at any point is the intersection of the latitude plane and the horizon plane at that point. Thus to 
obtain the east-west line we can take the cross product of the normals n^, to these two planes. 


From Figure 1, we can choose = (sin a, cos a, 0), and from Figure 2 we can choose = (0, BC, OC), so : 


n^xn^ = 


i j k 

sin a cos a 0 
0 BC OC 


parallel to 


= (OC cos a,—OC sin Of, i?C sin Of), 

, .BC 

(cos a, — sm a, -r-^ sm a) 

C/Cx 


From Figure 1 the easterly direction at point B must have a positive x-coordinate and a negative y-coordinate and 
thus we can choose the following as east and west (non-unit) vectors at B : 

( .BC. . BC. 

e = (cos a, — sma, sin a), w = (— cos a, sin a, — 7^77 sin a) 

C/O L/C/ 

Then taking the mirror image of these in the xy-plane, which reverses the 2 -component, with east mapping onto west 
and vice-versa, east/west vectors at B' are : 

, t .BC. . BC. 

e = (—cosQf,sma,sma), w = (cosa, — sma, ———sma) 


OC 


OC 


From the assumptions 3 and 5 the rays from the sun striking every point of the planet are parallel and come from the 
s = — i direction, and it is clear geometrically that the direction of s in the horizon plane at B' is south of the easterly 
e' at B'. Thus the desired magnitude of angle 9 of sunrise satisfies : 


e' • s = je'l jsj cos9, 
cosa 


.'. cos 9 = 


But 


/,2 . BC^ . ^ 

H = 1 + ^ sm a 


and so we need to express BC/OC in terms of a and S. From Figure 2, OC = RP^ — BC^, and BC equals the 
distance OB from Figure 1. From the triangle AAOB of Figure 1, we have : 


sin(90 — a — pi) sin(90 -I- a) . cos(a -I- pi) cos a 
— , le. — 


OB 


R 


OB 


R 
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Figure 1: Winter Solstice 



Thus BC of Figure 2 satisfies^ : 


BC = 


i?cos(a + f3) 
cos a 


Rs'mS 

COSO ’ 


and so : 


OC 


BC 

OC 


|e? = l + 


BC^ 

OC^ 


sin^ a 



sin (5 


^ V cos^ a — sin^ 6 

cos a 


\J cos^ a — sin^ d 

(cos^ a — sin^ <5) + sin^ d sin^ a cos^ a — sin^ <5 cos^ a 
cos^ a — sin^ S cos^ a — sin^ S 


cos^ a cos^ 6 
cos^ a — sin^ <5 


=> cosQ 
=> sin^ Q 
^ sind 


\J cos^ a — sin^ <5 


cos <5 


1 - 


cos a — sin 


cos^ d 



sin a 
cos (5 


1 — cos^ a sin^ a 

cos^ d cos^ (5 


^The special case of a = 90° can be checked separately using a plan view. 
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And thus as 9 above was an angle south of east the required angle north of east is given by equation (6). 


Figure 2: Winter Solstice - View Towards Sun 



V 


2.1.2 Southern Hemisphere 

The case S G [a — 90,0) in the southern hemisphere is shown in Figures 3 and 4. 

From the triangle AAOB of Figure 3, we have : 

sin(a + 13 — 90) sin(90 — a) 

OB ^ R 

— cos(a + /3) cos a 

OB ^ R 

Thus the y-coordinate BC of Figure 4, which is negative of OB of Figure 3, satisfies : 

i?cos(a + /3) RsinS 

BG = -= -, 

cos a cos a 

as before — thus we obtain the same east-west line as above and since the easterly direction at point B still must 
have a positive ^-coordinate and a negative y-coordinate the same expressions for e' and w', and hence 0, are obtained 
as in §2.1.1, thus again giving equation (6). 
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Figure 3: Winter Solstice - Case a + /3 > 90' 



Figure 4: Winter Solstice - View Towards Sun - Case a + j5 > 90° 
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2.2 Summer Solstice 


Here ip = 180° so we require to show 9 is minus what it is in the winter solstice case. The situation is as Figure 1 with 
the sun at the opposite side so s = i, and then it is clear a sunset 6 south of west becomes a sunrise 9 north of east, and 
a sunrise 9 south of east becomes a sunset 9 north of west, as required. 


2.3 General Day of Year 

Consider a day in the first quadrant Q1 as shown in Figure 5. The xyz-axes are fixed to the planet, but have constant 
direction as the planet orbits (note this is a different xyz-axis system than was used in Figures 1-4 of §2.1). The vertical 
planes through the solstice and equinox axes intersect at right angles in the vertical line of the orbital axis passing through 
the sun [18]. By rotating the NS axis by some angle p about vector b into plane F the effective rotation axis N'S' is 
produced, making effective axial tilt angle ai with the vertical. This is equivalent to simply viewing the sun and planet 
from a different angle, by rotating ourselves backwards by p from plane F. After that rotation the situation is then 
exactly the same as the winter solstice of Figure 1, except with the axial tilt angle ai instead of a. Thus the required 
angle 9 is given by equation (6) with a replaced by ai. Thus to obtain formula (1) we need to show sinoi = sin a cos ^/>. 


Let a = sinai + cosak be the NS axis unit vector, and let a' be the rotated N'S' axis unit vector. Then the required 
angle ai satisfies : 


cos oi = a' • k 

Wrt right hand rule the rotation about b = cos ip i + sini/; j is by —p, thus from the Rodrigues Rotation Formula 
(Appendix D) : 

a' = (cos p)a+ (sin p) (a x b) + (1 — cos p) (b • a) b 
Viewing along the b axis the rotation angle p satishes : 


a^^ • k = I a^ I cos p 


where a^ is the component of a perpendicular to b. Since a • b = sin a cos i/', 

a^^ = a — (a • b)b 

= a — (sinQ;cos'0)b 

= sin a sin^ ip i — sin a sin ip cos ip i + cos a k 

Thus 


I a^^ P = sin^ a sin^ ip + sin^ a sin^ ip cos^ ip + cos^ a 
= sin^ a sin^ ip + cos^ a 
= 1 — sin^ a cos^ ip 


Since a^ • k = cos a we then have (noting the case of the denominator zero (=> ■0 = 0) has already been covered) : 

COSO 

cos p = = 

a/I — sin^ a cos^ ip 

, . 2 1 cos^ a 

and, sm p = I- 2 -^— 

I — sin a cos^ ip 


Using 


1 — sin^ a cos^ ip — cos^ a 
1 — sin^ a cos^ ip 
sin^ a sin^ ip 
1 — sin^ a cos^ ip 


a X b 


i j k 

sin a 0 cos a 

cos Ip sin Ip 0 

— cos a sin ip i + cos a cos ip j + sin a sin ip k 


and the k-components of the above Rodrigues formula for a' we obtain 


a' • k = cos p cos a + sin p sin a sin ip 



Figure 5: General Day Of Year 


ORBITAL AXIS 



Then, using the above expressions for cos p and sin p, 

9 / / . x 9 (cos^ a + sin^ a sin^ ibY 

cos^ ai = (a' • k)^ = ^^ ' 

1 — Sin a cos^ yj 

(I — sin^ ocos^ V')^ , .9 9 , 

= - 7-5 ---= 1 — sin a cos 1/1 

1 — sin a cos^ ijj 

sin^ ai = sin^ a cos^'i/' 

sinai = sin a cos "0 

as required. 

For quadrant QA, the situation at day of year angle (360 — 1 /’)° is exactly the same as for Q1 at day of year angle i/^, 
with the same effective axial tilt ai, just on the opposite side of the solstice axis. Thus 9 must be the same, and then as 
cos(360° — tp) = costp, the formula (1) gives the correct result for QA. 
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For quadrant QS, the situation at day of year angle (180 + ip)° is exactly the same as for Q1 at day of year angle ^ 
except with the sun in the exact opposite direction, and thus as with the case of deriving the summer solstice case from 
the winter solstice case above the sunrise direction 6 just changes sign. Since cos(180° + ip) = — cosip, formula (1) then 
gives the correct sunrise angle for Q3. 

Finally for Q2, the situation at day of year angle (180 —■0)° is exactly the same as for Q3 at day of year angle (180+^)°, 
just on the opposite side of the solstice axis. Thus 9 must be the same, and then as cos(180° — ip) = cos(180° + ip), 
formula (1) gives the correct sunrise angle for Q2. 


2.4 Properties of Curve 

For a given a and S the equation (1) in radians has form f{x) = — sin“^(fccosa;), for x £ [0, 27r], where k 
and k £ (0,1) for non-polar latitudes (ie |5| < — a). Then : 


f'(^) 

and f"(x) 


k sin a: 

^/{l — k"^ cos^ x) 
fc(l — k^) cosx 
(1 — cos^ 


Thus the curve is ‘smooth’ at non-polar latitudes, with the following ‘bell shape’ gradient : 


sin a/ cos (5, 


X 

0 


•K12 


TT 


37r/2 


27r 


0 

+t 

k 


0 


-k 

-t 

0 


At a polar circle (ie 6 = ±(57r — a)), k = 1 and so in Q1 {x G [0, f{x) = — sin“^(cosa;) = x — ^tt, ie a linear 
function. Then by reflecting into Q2, Q3, and Q4 as above the overall function is a triangle. In the plots in Appendix A.l 
and in Figure 6 below it can be seen how this triangle is approached at the higher latitudes such as Reykjavik at 64.15°. 
Beyond the polar circles, fc > 1 and f{x) is undefined at certain times of the year when there is perpetual day or night, 
namely where | cosx| > 1/k. 


Figure 6: Sunrise Direction Curves 


fix) 



Orbital Position ip 
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3 Solar Declination Formula 


The solar declination angle [19] A is the angle made by the direction of the sun from the planet’s center with the planet’s 
equatorial plane, one of the two coordinates of a celestial body within the equatorial coordinate system [20]. From 
Figure 1 it is clear that A = — a on the winter solstice. And throughout Q1 since the situation of the planet is as Figure 1 
with a replaced with the effective axial tilt ai, we must have A = —ai. Thus in Ql, A € [—a, 0] and : 

sin A = —sin a cos f/' (7) 

In Qd, at angle (360 — 'tp)°, by symmetry A is the same as at angle xp in Ql, and as cos(360 — xp) = cosxp, formula (7) 
thus gives the correct result for QA when xp is substituted by (360 — xp)°. 

In Qi at (180 + xp)° the sun is in the exact opposite direction in space wrt planet’s center as it was in Ql at xp, and 
the solar declination is —A G [0,a]. But from (7) : 

sin(—A) = — sin a cos(I80 + xp) 

so that (7) gives the correct solar declination for Q3. 

In Q2, at (180 — xp)°, by symmetry solar declination is as Q3 at (180 + xp)°, ie. —A G [0, a]. But from (7) : 

sin(—A) = — sinacos(180 — xp), 
thus (7) gives the correct solar declination for Q2. 

Thus throughout the orbit : 

A = — arcsin(sina cosxp), A G [—a, a], (8) 

and the sunrise direction formula can be written as : 

0 = arcsin (, 9 G [—90°, 90°] (9) 

\coso / 

The formulae in the sections below are expressed in terms of A and equations (7) and (8) relate A to the day of the year 
via xp. The relationship of xp to day of the year depends on the orbital model used, for example the two orbital models 
considered in §1. 

4 Sunrise Equation and Daylight Duration Formula 

Consider a day in Ql as depicted in Figures 1 and 2 with a replaced with the effective axial tilt ai. The latitude circle is 
shown face on in Figure 7. A is the point of solar noon on the latitude circle, B is the point of sunset, and B' is the point 
of sunrise. The angle cr is the sunrise (and sunset) solar hour angle [21] which is a measure of the time between the sunrise 
(or sunset) and solar noon [22]. The time of solar noon is uniquely defined for any day of the year and for any arbitrary 
location L on the planet except for the poles which are stationary throughout the day^ (in the simplihed model). It is 
defined as the time when the unique half plane containing the location L and hinged on the NS axis intersects the sun — 
solar midnight is the time when this half plane is in the exact opposite direction from that. The general angular position 
of the half plane through L is a measure of time at L throughout the day as the planet rotates at a constant angular 
velocity about the NS axis. The solar hour angle r of T at a given time is the angle of Us half plane wrt the solar noon 
position, and is an angle in the range [0°,360°), with the +ve angle direction defined by the RH rule wrt axis NS. The 
solar hour angle is 0° at solar noon and 180° at solar midnight, and is negative at sunrise and positive at sunset. If the 
sun were placed in the exact opposite direction then the solar hour angle wrt the new sun position would be the previous 
solar hour angle plus 180° (this is used in §7 below). 

The latitude —d in the southern hemisphere is obtained from the latitude +<! by rotating Figure I about O by 180° 
so in Figure 7 the day/night arc lengths in (ii) are just the reverse of those in (i). Considering case (i) of northern 
hemisphere the angle cr is in the range [0, 90°] depending on ai and <5, eg it is 90° for all 5 at the spring equinox (12 hours 
of day and 12 hours of night) and is 0° at the artic circle on the winter solstice where solar noon, sunrise, and sunset all 
coincide. If Rl is radius of latitude circle then from Figure 1 : 

= sin(ai + /3) = cos 5 
= R cos 5 


Rl 

Rl 


®To define solar noon at a pole a meridian longitude would have to be selected for the pole (see §8 


re cardinal directions at the poles). 
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Figure 7: Latitude Circle Hour Angles for Day/Night 



(i) Latitude <5 > 0, Northern Hemisphere (ii) Latitude —5 < 0, Southern Hemisphere 


and from Figure 7 we have sintr = But from Figure 2, = OC, and from §2.1.1 : 


OC 

R 

j cos^ ai — sin^ S 

cosai ' 

. sin a 

OC 

\J cos^ a\ — sin^ <5 

Rcos6 

cos Qfi COS (5 

cos^ a 

cos^ ai cos^ S — cos^ ai + sin^ S 


cos-^ ai cos^ S 

, COS a 

= tan ai tan 5 

= — tan A tan 6 

(T 

= arccos(- 

- tan A tan 5) 


• 2 c- • 2 

sm 0 sm «! 2 c 2 

- - - — = tan 0 tan ai 

cos^ ax COS"' 0 


( 10 ) 


The identity (10) also applies in case (ii) of the southern hemisphere, at latitude —5, since from Figure 7(ii) the 
required sunrise solar hour angle is 180° — cr, and arccos as a function from [—1,1] to [0,180°] satisfies the condition 
arccos(—a;) = 180° — arccos(a;). The equation (10) is called the ‘Sunrise Equation’ [23], [24]. The argument to the arccos 
will go out of range for latitudes beyond the effective polar circles (ie the polar circles corresponding to the effective axial 
tilt), for at these latitudes and times there is no sunrise nor sunset. 


We have shown (10) holds for any day in Ql, using Figures 1 and 2 with an effective axial tilt ai G [0,a]. For a day 
in Q4 at angle (360 — ^/>), by symmetry the situation of the planet is the same as for a day in Q1 at angle ip, thus a is 
the same. Then since cos(360 — ip)° = cosip and from (8) the solar declination A is the same, equation (10) produces the 
correct sunrise solar hour angle a for QA. 


For a day in Q3 at angle (180 + ip)° and with solar declination —A, the geometric situation of the planet is identical 
to a day in Ql at angle ip with solar declination A (< 0), except with the sun’s light coming from exactly the opposite 
direction from that in Figure 1. For the day in Ql, equation (10) gives the value of a with the sun at its normal side 
of Figure 1, and then the required sunrise solar hour angle for the sun at the opposite side is (180 — cr)°, since the day 
and night arcs on the latitude circles will be reversed. Thus the required sunrise solar hour angle for the day in Q3 is 
(180 — cr)°. But from (10), for any latitude S, we have : 

180° — a = arccos(tan(—A) tan 5) 
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and so (10) gives the correct sunrise solar hour angle for the day in Q3 of solar declination -A. 

For a day in <32 at angle (180 — by symmetry the situation of the planet is the same as for a day in Q3 at angle 

(180 + and so a is the same. Then since cos(180 — 'tp)° = cos(180 + ^/>)° and the solar declination A from (8) is the 
same, the equation (10) produces the correct sunrise solar hour angle a for Q2. 

Thus (10) at any point in the orbit gives the required sunrise solar hour angle a G [0,180°]. Writing tan A in terms of 
the orbital angle 'll), (8) and (10) give : 


, , sinA sina cos-ip 

tan A = ±—- 2—7 =- 77 - 2 -7—7 (noting declination is -ve m Ql and CJA) 

•^(l —sin A) V(l —sin acos^ip) 

f sina cos if) 

and a = arccos I tan o ■ — - 2 - 7—7 

\ V(l ~ sin a cos^ i/>) 

To convert the total solar hour angle 2 ct of daylight to the daylight duration D we simply multiply the rotational 
period T of the planet by the fraction 2(j/360° = (t/180° : 



D = T 


a 

180° 


( 12 ) 


Although we are using the model of a stationary day, we could use as T the mean solar day [25], [26], of 24 hours 
duration in the case of Earth, rather than the sidereal day [27] length of 23.93447 hours, to compensate for the prograde 
motion of the Earth in its orbit during the day, which means it takes slightly longer than a sidereal day for the sun to 
return to the local meridian position at solar hour angle 0°. A planet with a retrograde rotation would have a mean solar 
day shorter than the sidereal day. 


In the case of Earth two sources of inaccuracy in equation (12) cause underestimation of the daylight duration D : 

• The non-zero angular diameter [14] of the solar disc, approximately 0.5°, and the convention of measuring sunrise 
and sunset wrt the upper edge of the solar disc crossing the horizon means official sunrise time comes a few minutes 
before and official sunset time a few minutes after the times predicted by the simple geometric model. 

• The effect of atmospheric refraction [28] of light from an object in space arriving at an observer on Earth causes 
the altitude of the object appear slightly higher than it really is. This is because entering the slower medium of 
the atmosphere from the vacuum causes the light ray to bend towards the normal between the two media, ie the 
angle of incidence reduces, and this implies a downwards bending, thus the eye sees the light source as higher up, 
similarly to a stick under water appearing higher up. This means the Sun can be seen shortly before sunrise and 
shortly after sunset, whilst its upper edge is still geometrically below the horizon, again lengthening the day. 

[3] provides details on adjustments which can be applied to (12) to compensate for these two factors. 


Note: the two sunrise angles a and 6 should not be confused — ct is the angle of the sunrise wrt solar noon, whilst 6 
is the angle of the sunrise wrt due east on the horizon — the former is an angle in the latitude plane and the latter 
is an angle in the horizon plane. Both are functions of time of year and latitude and from the equations (9) and (lO) 
cos (t/ sin 0 = — sin (5/cos A, which varies throughout the year. 


5 Clock Time Versus Solar Time 

In this section we relate the clock time CT (or civil time) to the solar time for the Earth^. The solar time is based on the 
daily passage of the Sun, with a day defined as the time between successive solar noons (§4). This is refered to as the 
apparent solar time (AST) or true solar time [29], and its day length varies with the time of the year, by about -20/-I-30 
seconds about a mean value of 24 hours. AST is the type of time measured by a sundial, as it is purely dependent on 
the position of the Sun which is recorded by the sundial’s shadow. The mean value of 24 hours forms the basis of the 
mean solar time (MST), and a mean solar day is of fixed duration 24 hours. As the successive daily differences accumu¬ 
late the time difference AST—MST varies within a range of about -I4/-I-16 minutes according to the ‘equation of time’ [30]. 

Clock time in the various time zones is based on the UTC time standard, which is within 0.9 seconds of another 
standard called UTl [31] which measures MST at longitude 0° (UTl replaces the older GMT standard). Thus unless 
greater than this level of accuracy is required UTC can be taken as the MST at longitude 0°, [32]. Each time zone is then 


"^We shall assume a 24 hour clock. 
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defined as UTC plus or minus an integral number of hours®, and has a defining meridian [32] located at longitude the 
corresponding integral multiple of 15°. The time zone time, ie the clock time CT, or civil time, shared by each location 
in the time zone, is defined as the MST at the longitude of the defining meridian for the zone. To determine the MST 
for any other location in the time zone we need to adjust by 4 minutes for every degree of longitude difference AL° from 
the longitude of the defining meridian for the zone, with the easterly direction being positive [32], ie. 

MST = CT + 4 minutes x AT° 

AL° 

MST = CT H—(in decimal hours) (13) 

1 0 

Daylight Saving Time (DST) is an adjustment to civil clock time of +1 hour for the summer months (eg from March 
to October in UK), but here we will define the civil clock time CT to be exclusive of any DST changes. 


When requiring the time of day in formulae (eg the solar hour angle parameters cr in the sunrise equation in §4 
and T in the altitude and azimuth formulae in §7 and §8) the AST is convenient, since geometrically it corresponds 
straightforwardly with the solar hour angle parameter, with 24 hours of AST being equivalent to a solar hour angle 360°. 
AST of 12 noon is solar noon, and AST of 12 midnight is solar midnight. An ‘hour’ of AST would be a 24*'' part of the 
solar day and would vary from day to day, within a range of approximately 3599 to 3601 atomic seconds. To obtain the 
time of day solar hour angle r from AST : 

(AST-12) = ^•24=§ 

.'. T = 15 (AST — 12), AST in decimal hours, t in degrees (14) 

To obtain r from the clock time CT we can use (13) and (14) together with the equation of time (EOT) which expresses 
AST in terms of MST, [30], [32], [33], and which is a continuous function of time throughout the year, repeating itself 
over a yearly cycle. We can approximate EOT as a constant for each day of the year d so that EOT(c?) equals a constant 
difference AST — MST throughout the day d. EOT((i) can thus produce AST from a known MST, and vice-versa. The 
difference EOT(c?) has two components : EC((i) the difference due to the eccentricity of the Earth’s orbit, and OB(d) 
the difference due to the obliquity of the ecliptic plane. EC(c?) is an approximate sine wave of period one year with zeros 
at perihelion and aphelion, which are respectively about 2 weeks after the winter and summer solstices (the exact day 
varying with the year), [34]. OB((i) is an approximate sine wave of period a half year with zeros at the equinoxes and 
solstices. Using respective amplitudes of 7.66 mins and 9.87 mins for these from [30], and taking a perihelion/aphelion 
offset from the winter/summer solstice of 14 days (the 2020 figure — other years will require a different value [34], [35]), 
and refering to the graphs in [30], [33], approximate formulae for the discretized EC((i) and OB((i) in decimal hours for 
the day which is at offset d from the winter solstice are : 

EC(d) = -0.1277 sin ( ■ 360° 

\ 365 

OB(d) = -0.1645 sin f• 720° 

\365 

and EOT((i) = EC(ci) -I- OB(c?) (in decimal hours) (15) 

Using the relation EOT((i) = AST — MST and equations (13) and (14) we then have : 

r = 15 (MST -h EOT((i) - 12) 

A r ° 

= 15 (CT + —— + EOT(d) - 12) 

15 

=> T = 15 CT — 180° -I- AL° + 15 EOT(d) (CT, EOT in decimal hours, r in degrees) (16) 

This is an approximate relation from clock time CT (ex-DST) to solar hour angle r for the day d, and is applicable to 
all time zones. Values for the approximate EOT((i) can be calculated using the above sine wave formulae or obtained from 
online calculators such as [36] and [37]. Note the approximation is based on a discretized form of the EOT — however it 
is sufficiently accurate to demonstrate the effect of the analemma and to generate approximate analemma graphs (§10). 


®In a few time zones, eg in India, a non-integral number is used. 
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6 Sunrise, Sunset, and Solar Noon Times 


The solar noon local clock time CTq for Earth is derived from equation (16) by putting r = 0 : 

AT° 

CTo = 12 - - EOT{d) 

1 0 


( 17 ) 


If (J G [0,180°] is the sunrise/sunset solar hour angle of §4 then in equation (16) r 
and r = cr at sunset clock time CTs- Then from (16) : 


CTfl = CTo-^ 

15 

CTs = CTo + ^ 

15 


—a at sunrise clock time CTr, 

(18) 

(19) 


with cr given by equation (10) or (11). 


Example Consider the location Madrid on the date 15**' May 2019. Longitude is —3.72° and latitude 6 = 40.42°. The 
time zone is UTC+1 (ex-DST) with defining meridian longitude at +15°. Thus AL° = —18.72°. Day offset from winter 
solstice 21’'* Dec 2018 is d = 10 + 31 + 28 + 31 + 30 + 15 = 145. Assuming a circular orbit/constant speed model to 
calculate an approximation to we have ip = (d/365) • 360° = 143.0°. Using EOT formula (15), EOT(d) = 0.05917 (the 
above EOT calculators [36] and [37] give values within about 3% of this). Then from (17) : 


CTo 
Actual CTq 


1 R 72 

12 + —-- 0.05917 = 13.19 hours = 13:11:24 (ex-DST) 

lo 

14:11:00 (DST) = 13:11:00 (ex-DST) 


and from equation (11) : 


^CT,^ 
and CTs 
Actual CTs: 
Actual CTs 


arccos 


/tan40.42 • sin23.44 • cos 143.0\ 
V C(1 - sin^ 23.44 • cos2 143.0) ) 


106.6° 

13.19- 106.6/15 = 6.083 hours = 06:04:59 (ex-DST) 
13.19 + 106.6/15 = 20.30 hours = 20:18:00 (ex-DST) 
06:58:00 (DST) = 05:58:00 (ex-DST) 

21:23:00 (DST) = 20:23:00 (ex-DST) 


where the actual sunrise, sunset, and solar noon times are obtained from [38]. The actual sunrise comes about 7 minutes 
before the time predicted by the model, and the actual sunset comes about 5 minutes after the time predicted by the 
model — this correlates with the two sources of discrepancy described in §4. The sunrise/sunset times given by [6] take 
into account the effect of refraction [39]. The solar noon prediction by contrast was more accurate and is not impacted 
by the above two factors. 


7 Solar Altitude Formula 

The solar altitude angle /i [40] is the angle the line from planet to sun makes with an observer’s local horizon plane. It 
is well-defined at all points on the planet at all times and is in the range [—90°, 90°], with the sun visible for /x > 0 and 
invisible for /x < 0, and with /x = 0 at sunrise and sunset. The formula derived here gives yx as a function of latitude 6, 
the day of the year (via A, as described in §3), and r the time of day expressed as the solar hour angle, §4 and [21], wrt 
solar noon [22]. Solar noon time can be looked up for any location and day in [6] or calculated using (17). r is negative 
before solar noon and positive after. 

Assume we are at a day in Q1 with effective axial tilt ai = —A, where A is the solar declination (§3) — ie. in 
Figure 1, a is replaced by ai — and consider a local horizontal coordinate system [41] at the point A in Figure 1 at 
latitude S with orthogonal right-handed basis vectors u (east), v (north), w (upwardly vertical) given (in terms of the 
a;j/ 2 :-coordinate system of Figure I) by u = (0, 0,1), v = (cos /3, sin/3, 0), w = (— sin/3, cos /3, 0). Note A could be any point 
on the half circumference of the planet including the poles N, S and /3 G [—a, 180° — a]. The other half circumference can 
be obtained by rotating about AS' by 180° — this ensures v always points north (see §8 re cardinal directions at the poles). 

The corresponding horizontal coordinate system basis vectors u' (east), v' (north), w' (upwardly vertical) for the 
point on latitude 5 of Figure 1 which is at solar hour angle t wrt solar noon can be obtained by rotating u, v, w by angle 
T about the polar axis NS according to the RH rule. Then at this angle t, if s = — i is the direction from the planet to 
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the Sun, making altitude angle ^ € [—90°,90°] with the u'v'-plane, then : 

^ > 0 s • w' = cos(90° — fi) = sin^ 

and fi <0 s • w' = cos(90 + |a^|) = — sin \^\ = sin ^ 

V/x G [—90°, 90°], sin^ = s • w' = —i • w' (20) 

Since in Figure 1 the NS rotation axis a = (sinai, cosai, 0), the Rodrigues Rotation Formula (Appendix D) gives : 
w' = (cost)w + (sinr)(a X w) + (1 — cosr)(a • w)a 

i j k 

axw = sinai cosai 0 = (0,0, sinai cos/3 + cosai sin/3) 

— sin /3 cos /3 0 

= (0, 0, sin(Q;i +/3)) = (0, 0, cos(3), and a • w = cos(q;i +/3) = sinJ 

w' = cosr(—sin/3, cos/3,0) + sinT(0,0, cosJ) + (1 — cost) sin5(sinQ;i, cosai, 0) 
sin fj, = — w' • i = cos r sin/3 — (1 — cos r) sin 6 sin ai 

Since 5 = 90 — (ai + /3), sin/3 = cos(ai + <3) 

sin/x = cosr(cosai cos5 — sinai sin 5) — (1 — cost) sin5sinai 
sin/x = COST cosai cos<3 —sin(3 sinai 

at any point in Ql, sin/x = cost cos <3 cos A + sin 5 sin A, /x G [—90°, 90°] (21) 

Note the above derivation is valid for any latitude on the planet, including the poles which are stationary throughout 
the day with <3 = ±90° and the formula (21) reducing to /x = A (for N) and /x = —A (for S) which is clear intuitively as 
the horizon plane at a pole is parallel to the equatorial plane. 

The above derivation was based on Figure 1 which covers case Q1 using the effective axial tilt ai in place of a in the 
diagram. In Q4 at (360 — ip)°, the situation of the planet wrt the sun is as Q1 at angle x// so /x must be the same. Since 

A given by (8) is also the same, the angle /x G [—90°, 90°] given by (21) is then correct for Q4, for any latitude S and time r. 

In Q3 the situation of the planet at angle (180 ± 7^)° is as Q1 at angle ip, except with the sun at the exact opposite 
side. Considering this symmetric Q1 position, with the sun at its usual side, as depicted in Figure 5, the solar hour angle 
is T± 180, where t is the solar hour angle for the Q3 situation (consider the half planes which define the solar hour angles 
and that the solar hour angle ±ve direction is always given from the NS axis via the RH rule — thus to go from a sun 
at one side to the complete opposite side add 180° to the solar hour angle), and the solar declination is —A where A is 
the solar declination for the Q3 position. Then, writing equation (21) for this Q1 point we have : 

sin/x = cos(t ± 180) cos(3 cos(—A) ± sinJ sin(—A), (22) 

and hence sin(—/x) = cos t cos 5 cos A ± sin (3 sin A (23) 

where /x is the solar altitude for the Q1 point, and —/x is the required solar altitude for the point in Q3. Thus equation (21) 
gives the correct solar altitude for the Q3 point. 

In Q2 the situation of the planet at angle (180 — ipy is as Q3 at angle (180 ± x/>)° so /x must be the same. Since 
cos(180 — Tpy = cos(180 ± x/>)° and A given by (8) is the same, the angle /x G [—90°, 90°] given by (21) is then correct for 
Q2, for any latitude 6 and time t. 

Thus equation (21) gives the solar altitude /x G [—90°,90°] at all latitudes <3 (including the poles), at all times of the 
year, and all times t of the solar day. The solar declination parameter A depends on the day of the year as described in 
§3. Equation (21) also gives an alternate route to the Sunrise Equation (10) by setting altitude /x = 0 in (21) and then 
solving for solar hour angle t. 


7.1 Sun Directly Overhead 

From equation (21), since the rhs is maximal at t = 0 the sun can only be directly overhead when t = 0, ie at solar noon. 
The condition for that to happen for a day is then : 

sin 90° = cos (3 cos A ± sin (3 sin A 

^ 1 = cos((3 — A) 

.-. as (3-A G (-180°,180°), 5-A = 0 

or sin(3 = sinA (as 5, A G [—90°, 90°]) 
from (7) sin(3 = —sinacosx// (24) 
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Outside the tropics, where |(5| > a, (24) has no solution in •0. On the Tropic of Capricorn at 5 = —a, there is a single 
solution 0 = 0, ie on the winter solstice. On the Tropic of Cancer at (5 = a, there is a single solution 0 = 180°, ie the 
summer solstice. Inside the tropical region : 


|(5| < Of 0 < 


sin S 
sin a 


< 1 


.'. there are two distinct solutions to (24) in 0 G [0, 360°), ie on two days of the year the sun is directly overhead at 
solar noon. In the case of the equator, these solutions are ip = 90° and ip = 270°, ie the two equinoxes. 


8 Solar Azimuth Formula 

Let the azimuthal angle in the horizon plane be 0 G [0°,360°) using the convention of measuring (p clockwise from due 
north [42]. Continuing from the main Q1 case in §7 above, let p be the component of s parallel to the horizon plane, 

ie parallel to the plane of u' and v'. Assume p 0, ie |/i| < 90°, otherwise <p is undefined. Then as (s • w')w' is the 

perpendicular component of s : 

p = s — (s • w')w' 

= —i + (i • w')w' 

= —i —w'sin/i , (from (20)) 

.’. IpP = P • P = 1 + sin^/r + 2(sin^)w'• i 

.’. jpj = cosfi (since from (20) sin/r = —i • w') (25) 

Now as (p is defined clockwise from north (ie from v') in the u'v'-plane, <p satisfies : 

jpj sin 0 = p-u', 

and jpj COS0 = p • v', 

.'. jpj sin 0 = (— i — w' sin/i) • u' = — u' • i (since basis vectors u', v', w' are orthogonal), (26) 

and jpj cos 0 = (— i — w' sin/r) • v'= — v'• i (27) 

But from the Rodrigues Rotation Formula : 

u' = (cost)u + (sinr)(a X u) + (1 — cosr)(a • u)a 

i j k 

axu = sinoi cosoi 0 = (cosai, — sinoi, 0), and a • u = 0 

0 0 1 

.'. u' = cosr(0, 0,1) + sinT(cosai, — sinai, 0) 

.'. u'• i = sinr cosai = sinr cos A 

and 

v' = (cost)v + (sinT)(a X v) + (1 — cosT)(a • v)a 

i j k 

axv = sinoi cosoi 0 = (0,0,sinai sin/3 — cosoi cos/3) 
cos (3 sin /3 0 

= (0, 0, — cos(ai +/3)) = (0,0, — sin0), and a • v = sin(ai +/3) = cos 5 

.'. v' = cos t(cos/3, sin/3,0) + sin t(0, 0, — sin0) + (1 — cost) cos5(sinai, cosoi, 0) 

.'. v'• i = COST cos/3 + (1 — cost) cos0 sin«! 

Since 0 = 90 — (oi + /3), cos/3 = sin(ai + 0) 

v'• i = COST sin(ai + 0) + (1 — cost) cos 0 sin «! 

= cos T(sin «! cos 0 + cos «! sin 0) + (cos t — 1) cos 0 sin ai 
= cos T sin 0 cos ot\ + cos 0 sin a\ = cos t sin 0 cos A — cos 0 sin A 

Thus from (25), (26), (27), (28), and (29) 

, — sinT cos A 

sm 0 = - 

COS/i 

, , cos 0 sin A — cos t sin 0 cos A 

and cos 0 = - 

COS/i 


(29) 

(30) 

(31) 
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Equation (31) can also be written without a r term : 


COSfj) 


(1 — sin^ 6) sin A — cos r sin 6 cos A cos 6 
cos /r cos S 

sin A — sin (5(cos r cos 6 cos A + sin <5 sin A) 
cos /r cos S 
sin A — sin S sin /x 


cos fi cos 6 


(from (21)) 


(32) 


Then equations (30) and (31) (or equations (30) and (32)) uniquely define (j) G [0°,360°), for the case of Ql. For the 
case of Q4, at an angle (360 — tp)° the situation of the planet is as at angle ip in Ql and thus the same (p applies. Since 
A given by (8) is the same at these two angles it follows that the rhs’s of (30) and (31) produce the same results, and 
hence give the correct (p for QA. 


For Q3 at angle (180 + 'tp)°, the same process of reflection as in §7 into Ql at angle ip can be performed, with the 
result that (30) and (31) produce negative of the Ql values for the sine and cosine of (p. Thus the required angle for Q3, 
which is {(p + 180°) mod 360°, is produced by (30) and (31). (Note the cos/x term does not change sign in the reflection). 


In Q2 the situation of the planet at angle (180 — ip)° is as Q3 at angle (180 + ip)° so (p must be the same. Since 
cos(180 — ip)° = cos(180 + ip)° and A given by (8) is the same, the angle (p G [0°, 360°) given by (30) and (31) is then 
correct for Q2, for any latitude S and time r. 


Thus equations (30) and (31) (or equations (30) and (32)) uniquely define the solar azimuth angle (p G [0°,360°) at 
all latitudes 5 at all times of the year, and all times t of the solar day, provided an azimuth is defined ie the sun is not at 
the zenith or nadir position. Since at the poles cardinal directions of north, south, east, and west are not well-defined® a 
meridian longitude (a half great circle) has to be selected, wrt which the cardinal directions can be defined. Then solar 
noon time is well-defined as the time when the sun crosses the selected meridian, even though on a stationary day in 
the orbit (in the simplified model) the sun’s altitude is constant (possibly negative) all throughout the day. North at a 
pole would then be defined as the continuation of northwards on the selected longitude line, which would be towards the 
north pole and away from the south pole. With this convention equations (30) and (31) produce the correct azimuth 
angle using the north clockwise convention, with r measured wrt to the selected meridian longitude : 


North pole ^ /x = A sin^ 
and cos (p 

<p 


— sinr 

— COST 

T -f 180° (mod 360°) 


which is the correct angle wrt due north as r is wrt due south, and 


South pole ^ /i = — A sin(() = — sinr 

and cos(p = cost 

(P = 360° - T (mod 360°) 
\e (p = —T (mod 360°) 


which is the correct angle wrt due north since at the south pole north points into the selected meridian half plane, and a 
positive rotation t from solar noon will produce a negative solar azimuth (p according to the north clockwise convention. 
The above process will work for any choice of meridian longitude for a pole. 


The azimuth cosine formula (32) provides an alternate route to the sunrise direction formula (9) since putting /x = 0 
in (32) and remembering 9 is measured anti-clockwise from due east we obtain : 

6 + (p = 90° sin 0 = cos (p = 

COSO 

This completely characterizes 9 since from §1 we know in all cases 9 G [—90°, 90°]. 


®At the north pole every direction points to the south and at the south pole every direction points to the north. 
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Example Continuing with the example of Madrid on 15*^ May 2019 in §6, compare the above altitude and azimuth 
formulae with the actual solar altitude of 50° and azimuth of 249° at DST 16:47 from [38]. We have : 

CT = 15:47 = 15.78 hours 

r = 15 • 15.78 - 180 - 18.72 + 15 • 0.05917 = 38.87°, .-. cost = 0.7786, sinr = 0.6275 
and sin A = — sin 23.44 • cos 143.0 = 0.3177, cos A = 0.9482 

sin/x = 0.7786 • (cos 40.42) • 0.9482 + (sin40.42) • 0.3177 = 0.7681 
/X = 50.2°, cos yx = 0.6403. 

And sincj) = -0.6275 • 0.9482/0.6403 = -0.9297 

cos(j) = (0.3178 - (sin 40.42) • 0.7681)/(0.6403 • cos 40.42) = -0.3697, from (32) 

^ (f) = 248.4° 

Thus the formulae predictions are accurate to within 0.6° when compared with [38]. 


9 Solar Noon Altitude 

The solar noon altitude [24] /kq for a location is by definition the solar altitude when the solar hour angle r for the location 
is 0°. From (21) this is the maximum value of /x throughout the day, ie the highest position of the sun. On a day of 
perpetual darkness this will be negative, eg as in Figure 1 at a point above the arctic circle on the winter solstice (where 
case 1 below applies). From Figure 1 with the effective axial tilt ai, covering cases Q1 and Q4, it is clear the azimuth (j) 
of the sun at solar noon for any day and latitude is either : 

1. southerly, ie ^ = 180°, for latitudes above the effective Tropic of Capricorn (at latitude —oi) 

2. northerly, ie </> = 0°, for latitudes below the effective Tropic of Capricorn 

3. undefined on the effective Tropic of Capricorn, as the sun is then directly overhead at solar noon (ie yxo = 90°) 

The effective Tropic of Capricorn is the southern tropic corresponding with the effective axial tilt ai. The effective 
tropics converge to the equator at the spring equinox where ai becomes zero. Note cases (1) and (2) do not correspond 
with the northern and southern hemispheres, but with the two sides of the effective tropic which changes throughout the 
year. For example in Figure 1 a southern latitude between —ai and 0° will see the noon sun in the southern half of the 
sky, as a northern latitude will. The dividing line is the effective tropic, which at the spring equinox becomes the equator. 
For Q2 and Q3 a similar set of three cases applies involving the effective Tropic of Cancer in the north at latitude ai. For 
practical purposes which of the above three cases applies would be determinable whenever the sun is visible at solar noon. 


Considering case 1, entering r = 0, </> = 180, /x = y:io in equations (21) and (31) we obtain : 


sm /io 
cos /io 
ie. cos(90 — /io) 
and sin(90 — /xo) 
(90 - /xo) -id-X) 


cos 6 cos A + sin (5 sin A 
sin S cos A — cos S sin A 
cos((5 — A) 
sin((f — A) 

integral multiple of 360 


But then, defining the solar zenith angle [40] vq at solar noon to be the complement of angle /xo : 


/io G [-90,90], (5 G (-90,90), A G [-a, a], a G [0,90) 
(90 -/Xo)- (5-A) G (-180,360) 

(90 - /io) - (5 - A) = 0 
vq = 6 — X 

Considering case 2, entering t = 0, </> = 0, /x = /xq in equations (21) and (31) we obtain : 


sin /io 
cos Ho 
ie. cos(90 — Ho) 
and sin(90 — ho) 


cos 6 cos A + sin 5 sin A 
cos 6 sin A — sin S cos A 
cos{6 — A) 

— sin((f — A) 


(33) 


(34) 


Thus in the xx/-plane on the unit circle the angles 90 — ho and 6 — X are reflections of one another in the x-axis. 
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90-fio = 360 - (,5 - A) (mod 360) 

270 + /io — + A = integral multiple of 360 

but from (33), - S + X G (-270, 270) 

270 + ^o-'5 + A G (0,540) 

270 + fXo-S + X = 360 

t'o = X — S (35) 

In case 3, entering r = 0, /io = 90 in equation (21) we obtain : 

sin 90 = cos 6 cos A + sin 5 sin A 

ie. 1 = cos{6 — A) 

S — X = integral multiple of 360 
from (33), S — X = 0 

^ X = 6 (36) 

The three solar noon equations (34), (35), and (36) give a relationship between latitude A and solar declination 
S depending on which of the above cases 1-3 applies. These equations will still apply even if fiQ < 0, ie a day of 
perpetual darkness, however in case 3 the one possibility not considered, ie /xq = —90, with the sun at the nadir 
position at solar noon, is not possible, as we would expect intuitively, but formally by (21) such a position would imply : 
sin(—90) = —1 = cos((l — A) => d — A — 180 = multiple of 360. But since from (33), S — X — 180 G (—360,0), this is 
impossible. The solar noon equations can allow the latitude in daytime to be determined from a table of solar declinations 
(which can be computed from (8)), by measuring the zenith angle of the sun at noon. 

10 Analemma 

One form of analemma [43] is a graph of the Sun’s position in the sky as seen from a fixed location on Earth, as altitude 
plotted against azimuth, at the same local clock time CT (ex-DST) each day for a period of a year. With equal horizontal 
and vertical scales this will have the same shape as a trace of the Sun’s actual positions over these times, for example as 
captured in analemma photographs. The irregular figure 8 shape of the analemma (Figure 8) is a result of the deviation 
throughout the year of the AST from the MST, as given by the equation of time (§5). In the present simplified model 
the analemma is given by equations (21), (30), and (31) together with equation (16) which for a fixed daily local clock 
time CT, produces a varying solar hour angle r for the day d, which is then entered in equations (21), (30), and (31) 
to produce the altitude /x and azimuth (p. The parameter A also varies daily, according to equation (8). If the solar day 
length were constant throughout the year so that MST = AST, and EOT{d) = 0 Vd, then t would be constant for clock 
time CT each day, and the analemma would be a simple curve without any loops because the solar declination parameter 
A retraces its values in [0,180°] in the return journey from 180° to 360° due to the costp term in equation (8). 

The Perl script calc-analemma.pl in Appendix C.3 calculates the analemma for a given location and time using 
the above formulae, producing a table of altitude and azimuth values over a year that can be directly input into the 
TIKZ/pgfplots package to produce a graph (eg by compiling the file analemma-graph.tex). The script uses the circular 
orbit/constant speed model to calculate an approximate orbital angle ip for each day of the year. To use the script specify 
the desired local clock time CT (ex-DST), the latitude S, and the difference AL° in degrees between the longitude of the 
location and the longitude of the defining meridian for the time zone^. 

Uncommenting a line in the script shows the effect of a constant solar day length, and other lines can be uncommented 
to show the separate effects of the eccentricity and obliquity components of the EOT. 

Note that an analemma in the tropics will become ill-defined towards solar noon because at such locations at certain 
times of the year the Sun is directly overhead at solar noon (§7.1) and the azimuth is then no longer well-defined. As we 
approach solar noon at such locations the azimuth equations (30) and (31) approach 0/0. Note because of equation (21), 
which implies /x is maximal at r = 0, the Sun directly overhead can only occur at solar noon (and likewise the nadir 
position only at solar midnight). 


^Note the longitude itself is not sufficient information as a single longitude can span more than one time zone, eg the UK and Spain span 
time zones UTC and UTC+1. 
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Since analemma photographs are rare the accuracy of the analemmas produced by the script are best checked by com¬ 
paring the overall accuracy of the altitude and azimuth formulae (21), (30), and (31) against real data. In the example 
of §8 the accuracy was 0.6°. 

In Figure 8 below two analemmas generated by the script are shown : 

(i) Athens, Greece : latitude 37.98° N, longitude = 23.73° E, UTC-b2, AL° = -6.27°, CT = 16:00 
Script command line : perl calc-analemma.pl 37.98 16 -6.27 > analemma-athens.dat 
Photo : http://www.perseus.gr/Astro-Solar-Analemma-140000.htm 

(ii) Kumagaya, Japan : latitude 36.15° N, longitude = 139.38° E, UTC-l-9, AL° = 4.38°, CT = 7:00 
Script command line : perl calc-analemma.pl 36.15 7 4.38 > analemma-kumagaya.dat 
Photo : https://earthsky.org/todays-image/todays-image-analemma-2013 


(i) Athens, Greece, 16:00 hours 


(ii) Kumagaya, Japan, 07:00 hours 




Azimuth (f> (°) 


Figure 8: Analemmas Generated by Formulae (16), (21), (30), and (31) 


For a planet other than the Earth, with an axial tilt a £ [0, 90°), provided the original assumptions 1-6 in §1 apply, 
the equations (21), (30), and (31) would still apply as they are in terms of solar time, but equation (16) for clock time, 
and the EOT, would be different and the mean solar day length different from 24 hours. Also if the orbit were not well 
approximated by a circle then a more complicated mapping from day of the year d to orbital angle tp would be required. 
With the required changes the script should then be able to generate analemmas for the planet. 
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Appendix A 


Comparison With Actual Sunrise Data 


A.l Sunrise Predictions 

The actual sunrise directions 9 for locations at 8 different latitudes (4 northern hemisphere, 3 southern hemisphere, 1 equa¬ 
torial) at 10 day intervals over the year 2018-2019 are shown in Table A.l below (source : https: //www. timeanddate. com, 

[6]). WWW .timeanddate . com gives the sunrise and sunset directions to the nearest degree. Where in the www. timeanddate. com 
figures the sunrise angle north of east differs from the sunset angle north of west on a day the average of these two is 
quoted in the table below. The day offset n is from the 2018 winter solstice date of 21/12/2018. 


Day Offset 
n 

Abu Dhabi 
24.45° 

Edinburgh 

55.95° 

Melbourne 

-37.81° 

Milan 

45.47° 

Quito 

-0.17° 

Reykjavik 

64.15° 

Rio 

-22.88° 

Stanley 

-51.69° 

0 

-25 

-44 

-31 

-34 

-23 

-62 

-26 

-41 

10 

-25 

-43 

-31 

-33 

-23 

-60 

-26 

-41 

20 

-24 

-40 

-29 

-31 

-22 

-56 

-24 

-38 

30 

-22 

-36 

-27 

-28 

-20 

-49 

-22 

-35 

40 

-19 

-31 

-23.5 

-25 

-18 

-42 

-19.5 

-30.5 

50 

-16 

-25.5 

-19.5 

-20 

-15 

-33.5 

-16 

-25 

60 

-12 

-19 

-15 

-15 

-11 

-25 

-12.5 

-19.5 

70 

-8 

-12.5 

-10.5 

-10 

-7.5 

-16 

-8.5 

-13 

80 

-4 

-5.5 

-5.5 

-4.5 

-4 

-7 

-4 

-7 

90 

0.5 

1.5 

0.5 

1 

0 

2.5 

0 

-0.5 

100 

5 

8.5 

4.5 

7 

4 

11.5 

4 

6 

no 

9 

15.5 

9 

12 

8 

20.5 

8 

12 

120 

13 

22 

14 

17.5 

12 

29.5 

12 

18 

130 

17 

28.5 

18 

22 

15 

38 

16 

23 

140 

20 

34 

22 

26.5 

18 

46.5 

19 

28 

150 

22 

39 

25 

30 

20 

54.5 

21 

32 

160 

24 

43 

27 

33 

22 

61.5 

23 

35.5 

170 

26 

46 

29 

35 

23 

67.5 

25 

38 

180 

26 

47 

29 

36 

23 

70 

25 

39 

190 

26 

46.5 

29 

35 

23 

69 

25 

38 

200 

25 

44.5 

28 

34 

22 

64.5 

24 

36.5 

210 

23.5 

41 

26 

31.5 

21 

57.5 

22 

34 

220 

21 

36.5 

23 

28 

19 

50 

20 

30 

230 

18 

31 

20 

24 

16 

41.5 

17 

25 

240 

15 

25 

16 

20 

13 

33.5 

14 

20 

250 

11 

19 

12 

15 

10 

24.5 

10 

14.5 

260 

7 

12 

7 

9.5 

6 

15.5 

6 

9 

270 

3 

5 

2.5 

4 

2 

7 

2 

2.5 

280 

-1 

-2 

-2.5 

-1.5 

-2 

-1 

-2 

-4 

290 

-5.5 

-8.5 

-7.5 

-7 

-5.5 

-11 

-6 

-10 

300 

-10 

-15.5 

-12 

-12 

-9 

-20 

-10.5 

-16 

310 

-14 

-22 

-17 

-17.5 

-13 

-28.5 

-14 

-22 

320 

-17 

-28 

-21 

-22 

-16 

-37 

-18 

-27.5 

330 

-20 

-33.5 

-24.5 

-26 

-19 

-45 

-21 

-32.5 

340 

-23 

-38 

-27.5 

-30 

-21 

-52 

-23 

-36.5 

350 

-24 

-41 

-30 

-32 

-22.5 

-58 

-25 

-39.5 

360 

-25 

-43 

-31 

-33 

-23 

-61 

-26 

-41 

365 

-25 

-44 

-31 

-34 

-23 

-62 

-26 

-41 


Table A.l: Actual Sunrise Directions 9 for 2018-2019 


The charts below compare the actual sunrise directions shown in red with the curve of sunrise directions predicted 
by the sunrise direction formula (with the currently accepted Earth axial tilt a = 23.44°) shown in blue. The circular 
orbit/constant speed model is used with the orbital angle in the formula calculated as -^(360°) for day offset d G [0, A^—1] 
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Degrees [9) North of East Degrees [9) North of East 


from the day of the winter solstice, where N = 365 is the number of days in the year. The average error over the 8 charts 
is 1.25°. When the formula’s own axial tilt estimate of 23.52° (Appendix B) is used the average error is 1.26°. 


Actual/Computed Sunrise in 
Edinburgh, 2018-19 



Day Offset (n) from Winter Solstice 


Actual/Computed Sunrise in 
Milan, 2018-19 



Day Offset (n) from Winter Solstice 
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Degrees {9) North of East Degrees [9) North of East 


Actual/Computed Sunrise in 
Reykjavik, 2018-19 



Day Offset (n) from Winter Solstice 


Actual/Computed Sunrise in 
Abu Dhabi, 2018-19 



Day Offset (n) from Winter Solstice 
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Degrees {9) North of East Degrees [9) North of East 


Actual/Computed Sunrise in 
Melbourne, 2018-19 



Day Offset (n) from Winter Solstice 


Actual/Computed Sunrise in 
Stanley, Falkland Islands, 2018-19 



Day Offset (n) from Winter Solstice 
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Degrees { 9 ) North of East Degrees [ 9 ) North of East 


Actual/Computed Sunrise in 
Rio De Janeiro, 2018-19 



Day Offset (n) from Winter Solstice 


Actual/Computed Sunrise in 
Quito, Ecuador, 2018-19 



Day Offset (n) from Winter Solstice 
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A.2 Perpetual Day/Night Predictions 

Using the same circular orbit/constant speed model as in Appendix A.l, consider the location of Jan Mayen, Norway at 
latitude 71° north. The argument to the arcsin goes out of range when either: 


cosi, > f2!l = _S££llL ^0.81845 = cos35" 

sma sm 23.44° 

cos S 

or cos'tp < -~ —0.81845. 

sma 


These correspond with ip ranges of [0, 35], [145, 215], [325, 360], which for the year 2018-2019 correspond with the dates: 
25/1/2019, 17/5/2019, 27/7/2019, 15/11/2019, for entering/exiting periods of perpetual day or night. The actual dates 
from timeanddate.com, [6] are: 


Actual date 22/1/2019 13/5/2019 1/8/2019 21/11/2019 

Estimated date 25/1/2019 17/5/2019 27/7/2019 15/11/2019 

Error (days) 3 4 5 6 

Average error 4.5 days 


For Cape Adare, Antartica at latitude 71° south, the estimated dates for 2018-2019 are as for Jan Mayen and the 
actual dates are: 

Actual date 31/1/2019 19/5/2019 26/7/2019 14/11/2019 

Estimated date 25/1/2019 17/5/2019 27/7/2019 15/11/2019 

Error (days) 6 2 11 

Average error 2.5 days 


For Longyearbyen, Svalbard, Norway at latitude 78° north, the arcsin goes out of range when cosip > 0.52267 ~ 
cos58.5° or cosip < —0.52267, ie in ip ranges [0, 58.5], [121.5, 238.5], [301.5, 360], which for 2018-2019 correspond with the 
dates: 18/2/2019, 23/4/2019, 20/8/2019, 23/10/2019. The actual dates are: 

Actual date 16/2/2019 19/4/2019 25/8/2019 27/10/2019 

Estimated date 18/2/2019 23/4/2019 20/8/2019 23/10/2019 
Error (days) 2 4 5 4 

Average error 3.75 days 

For McMurdo Station, Antartica at latitude 78° south, the estimated dates for 2018-2019 are as for Longyearbyen 
and the actual dates are : 

Actual date 20/2/2019 25/4/2019 19/8/2019 24/10/2019 

Estimated date 18/2/2019 23/4/2019 20/8/2019 23/10/2019 
Error (days) 2 2 11 

Average error 1.5 days 

The average errors overall is 3.1 days which though not as good as the sunrise predictions of §A.l or the axial tilt 
prediction of Appendix B are still reasonable approximations for the simplified model. 
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Appendix B 

Axial Tilt Estimates 


The axial tilt estimates of the sunrise direction formula obtained by applying equation (5) to the sunrise data of Table A.l, 
using the same circular orbit/constant speed model as in Appendix A.l, are shown in Table B.l below. The estimates 
can become abnormally inaccurate when observed values of sin 9 and cos 'ip are close to zero near the equinoxes (n ~ 90 
for spring and n ~ 270 for autumn), for then small absolute changes in these two terms can result in large % changes, 
causing large % change in sin a on Ihs of (5). The rhs of (5) is thus unstable near the equinoxes (and at the equinoxes it 
has form 0/0). The cos 5 term does not cause instability so long as we are not too close to the poles. A good direction 
approximation can have a large % error when we are near zero - eg approximating direction of 0.1° by 0.2° has 100% 
error but is still a good direction approximation. Thus the particularly bad 5 estimates near the spring equinox in the 
table which were negative are excluded in calculating the averages — though any other bad estimates are left in. Over 
the remaining 299 data points the overall average estimate for the Earth’s axial tilt is 23.52°, which is within 0.1° of the 
currently accepted value of 23.44°. 


Day Offset 

n 

Abu Dhabi 
24.45° 

Edinburgh 

55.95° 

Melbourne 

-37.81° 

Milan 

45.47° 

Quito 

-0.17° 

Reykjavik 

64.15° 

Rio 

-22.88° 

Stanley 

-51.69° 

Overall 

Average 

0 

22.63 

22.89 

24.01 

23.09 

23.00 

22.64 

23.82 

24.00 


10 

22.99 

22.80 

24.39 

22.81 

23.37 

22.54 

24.20 

24.38 


20 

23.16 

22.48 

24.01 

22.56 

23.45 

22.58 

23.46 

23.92 


30 

23.09 

22.24 

24.36 

22.25 

23.16 

22.24 

23.38 

24.14 


40 

22.57 

21.93 

24.08 

22.57 

23.59 

22.20 

23.47 

24.05 


50 

22.64 

21.70 

23.86 

21.59 

23.39 

21.66 

22.93 

23.70 


60 

21.68 

20.84 

23.52 

20.75 

21.86 

21.08 

22.90 

23.82 


70 

20.74 

19.80 

23.73 

19.90 

21.40 

19.63 

22.38 

22.95 


80 

19.27 

16.19 

23.17 

16.61 

21.25 

16.03 

19.51 

23.11 


90 

-21.67 

-42.94 

-18.69 

-34.67 

0.00 

-62.12 

0.00 

14.56 


100 

31.92 

33.47 

24.40 

34.72 

27.70 

35.40 

25.36 

25.58 


no 

26.68 

28.15 

22.93 

27.37 

26.03 

28.78 

23.84 

23.98 


120 

25.54 

26.21 

23.73 

26.36 

25.96 

26.88 

23.79 

23.79 


130 

25.48 

25.58 

23.24 

25.13 

24.73 

25.72 

24.23 

23.05 


140 

24.73 

24.88 

23.44 

24.87 

24.54 

25.15 

23.77 

23.02 


150 

23.73 

24.57 

23.20 

24.44 

23.80 

24.76 

22.93 

22.81 


160 

23.57 

24.36 

22.79 

24.36 

23.86 

24.45 

22.88 

22.88 


170 

24.11 

24.35 

23.08 

24.31 

23.58 

24.35 

23.49 

23.00 


180 

23.54 

24.20 

22.54 

24.37 

23.02 

24.21 

22.94 

22.98 


190 

23.73 

24.18 

22.72 

23.93 

23.20 

24.23 

23.12 

22.64 


200 

23.76 

24.26 

22.85 

24.25 

23.10 

24.34 

23.10 

22.71 


210 

24.07 

24.38 

22.90 

24.31 

23.74 

24.40 

22.82 

22.92 


220 

24.11 

24.64 

22.73 

24.34 

24.05 

24.72 

23.23 

22.83 


230 

24.29 

24.94 

23.27 

24.65 

23.77 

24.99 

23.19 

22.52 


240 

25.42 

25.54 

23.38 

25.91 

24.20 

26.01 

23.96 

22.73 


250 

25.91 

27.29 

24.41 

27.17 

25.90 

27.05 

23.73 

22.98 


260 

28.24 

29.77 

24.24 

29.58 

26.47 

29.80 

24.25 

24.43 


270 

47.61 

49.16 

32.29 

49.32 

32.75 

55.46 

29.90 

24.78 


280 

8.51 

10.48 

18.72 

9.84 

18.97 

4.06 

17.42 

23.75 


290 

18.42 

17.44 

21.93 

18.03 

20.31 

17.54 

20.41 

22.95 


300 

21.22 

20.04 

22.10 

19.51 

20.99 

19.97 

22.61 

23.04 


310 

22.14 

21.04 

23.29 

21.16 

22.64 

20.86 

22.42 

23.42 


320 

21.86 

21.58 

23.34 

21.57 

22.69 

21.54 

23.48 

23.61 


330 

22.20 

22.03 

23.43 

21.91 

23.27 

21.97 

23.62 

23.84 


340 

23.04 

22.29 

23.67 

22.69 

23.22 

22.21 

23.34 

23.94 


350 

22.52 

22.33 

24.11 

22.60 

23.32 

22.49 

23.75 

24.07 


360 

22.72 

22.54 

24.11 

22.54 

23.09 

22.51 

23.92 

24.09 


365 

22.63 

22.89 

24.01 

23.09 

23.00 

22.64 

23.82 

24.00 


Average*^ 

23.8 

23.88 

23.57 

23.9 

23.06 

23.98 

22.67 

23.29 

23.52 


Table B.l: Axial Tilt Estimates in Degrees (°) 


®excluding any negative estimates of axial tilt. 
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Appendix C 


Perl Scripts 

The scripts C.l and C.2 below use the circular orbit/constant speed model for Earth with stationary days evenly spaced 
around the circle and the orbital angle ip defined as ;^(360°) for day d £ [0, — 1], where d is the day offset from the 

day of the winter solstice, and N = 365 is the number of days in the year. To use the more accurate mapping from day 
of the year to ip described in [3] make the following changes to these scripts (the time of noon is used for each day, and 
each day duration is taken to be a mean solar day of 24 hours). A similar change can be made to script C.3. 

• replace the line my $winter_solstice = DateTime->new(day => 21, month => 12, year => 2018); with 
my $base_date = DateTime->new(day => 1, month => 1, year => 2013); 

• replace the line my $day_offset = $date->delta_days($winter_solstice)->in_units(’days’) ; with 
my $t = $date->delta_days($base_date)->in_units(’days’) + 0.5; 

• replace the line $psi = ($day_offset / 365) * pi2; ) with the formula of [3] as a function of $t, plus $pip2 (ie 
90°) 


C.l Calculate Sunrise Direction 


# Compare actual sunrise with sunrise computed using axial tilt of 23.44 degrees. Calculate average error. 

# Usage : perl calc-sunrise.pl <file 1> ... <file N> 

# Each input file line has format 

# <Day> <Month> <Year> <Actual Sunrise> <Elevation> 

# or 

# <Day> <Month> <Year> <Actual Sunrise> 

# where <Elevation> is the Pole Star (or southern equivalent) elevation above the horizon 

# (ie the ^latitude’ in the spherical model). 

# When elevation is omitted from a line the last specified value is used. 

# Lines containing only white space or beginning with are ignored. 

# Annotation lines begin with field ’aumot’ (case-insensitive) and are echoed to output. 

use strict; 
use warnings ; 

use 5.012_003; 

use DateTime ; 

use DateTime :: Duration; 

use Math::Trig; 

use Math::Trig qw(:pi); 

my $AXIAL_TILT = deg2rad(23.44); 

my $winter_solstice = DateTime->new( 
day => 21, 

month => 12, 

year => 2018 

); 


print ( "Actual Sunrise Computed Sunrise Abs Error\n"); 
printC - - - \n"); 

my $saved_elevation ; 
my $line_counter = 0; 
my $data_counter = 0; 
my $abs_total = 0; 

while (<>) { 

$line_counter++; 

my ©fields = split; 

my $fields_count = scalar ©fields; 

# ignore any line containing only white-space 
next if (! $fields_count ); 

# skip any line beginning with 

next if ( substr ($fields [0], 0, 1) eq ’#' ); 

# echo cinnotation line (eg a blaink line cam be inserted with this) 
if ( lc($f ields [0]) eq 'amnotO { 

shift ©fields; 
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> 


printf ("@fields\n") ; 
next ; 


# 4 cuid 5 are the 
if ($fields_count 
printf (STDERR 
exit 1; 

} 


only valid numbers of fields on a line 
!= 4 && $fields_count != 5) { 

"Invalid line : line number 7,d\n" , $line_counter ); 


my ($day, $month, $year, $actual_sunrise , $elevation) = Ofields; 


if ( ! defined ($elevation) ) { 

# only 4 fields were specified on the line - elevation not specified, so we use the last 

# specified elevation, exiting with an error if the latter does not exist 
if (! defined ($saved_elevation) ) { 

printf (STDERR "Input files do not specify elevation data on first data-containing line.\n") 
exit 1; 

} 

$elevation = $saved_elevation ; 

} else { 

# 5 fields specified on the line - save off the elevation 
$saved_elevation = $elevation; 

> 


my $date = DateTime->new( 
day => $day, 

month => $month, 

year => $year 


); 


# zero-based day offset from Winter Solstice 

my $day_offset = $date->delta_days($winter_solstice ) ->in_units ( Mays O ; 

# various angles 

my ($psi, $delta, $theta) ; 

# day of year angle 

$psi = ($day_offset / 365) * pi2; 

# calculate Pole Star Elevation in radiauis 
$delta = deg2rad($elevation) ; 

# sine of sunrise auigle 

my $sin_theta = -sin($AXIAL_TILT) * cos($psi) / cos($delta); 

# sunrise aingle 

$theta = rad2deg( asin($sin_theta) ); 

$dat a_ c ount er ++; 

# here Mbsolute^ means absolute error as opposed to relative error, not modulus 
my $abs_error = $theta - $actual_sunrise ; 

$abs_total += abs ($abs_error) ; 

printf ("7,5. If7,12sy,6.2fy,13sy,5.2f\n" , $actual_sunrise , $theta, $abs_error) ; 

} 

printf ("Average magnitude abs error = y,.2f degrees\n", $abs_total / $data_counter) ; 
printf ("7,d data points\n", $data_counter) ; 


Sample output : 


C:\>perl calc-sunrise.pl iedinburgh.txt 
Actual Sunrise Computed Sunrise Abs Error 


-43.0 -44.42 -1.42 

-40.0 -41.97 -1.97 

-36.0 -38.16 -2.16 

Edinburgh 

Average magnitude abs error = 2.14 degrees 
38 data points 
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C.2 Calculate Axial Tilt 


# Calculate axial tilt from actual sunrise data. 

# Usage : perl calc-axial-tilt.pl <file 1> ... <file N> 

# Each input file line has format 

# <Day> <Month> <Year> <Actual Sunrise> <Elevation> 

# or 

# <Day> <Month> <Year> <Actual Sunrise> 

# where <Elevation> is the Pole Star (or southern equivalent) elevation above the horizon 

# (ie the ^latitude’ in the spherical model). 

# When elevation is omitted from a line the last specified value is used. 

# Lines containing only white space or beginning with are ignored. 

# Annotation lines begin with field ’eumot’ (case-insensitive) and are echoed to output. 

use strict; 
use warnings ; 

use 5.012_003; 

use DateTime ; 

use DateTime :: Duration; 

use Math::Trig; 

use Math::Trig qw(:pi); 

my $winter_solstice = DateTime->new( 
day => 21, 

month => 12, 

year => 2018 

); 


printC'PSI cos(PSI) Actual Sunrise Axial Tilt\n"); 

printC - - - -\n"); 


my $saved_elevation ; 
my $line_counter = 0; 
my $data_counter = 0; 
my $total = 0; 

while (<>) { 

$line_counter++; 

my ©fields = split; 

my $fields_count = scalar ©fields; 

# ignore any line containing only white-space 
next if (! $fields_count ); 

# skip any line beginning with 

next if ( substr ($fields [0], 0, 1) eq ’#’ ); 

# echo annotation line (eg a bleuik line can be inserted with this) 
if ( lc($f ields [0]) eq ’annotO { 

shift ©fields; 
printf ("©fields\n") ; 
next ; 

} 


# 4 and 5 are the 
if ($fields_count 
printf (STDERR 
exit 1; 

} 


only valid numbers of fields on a line 
!= 4 && $fields_count != 5) { 

"Invalid line : line number 7,d\n" , $line_counter) ; 


my ($day, $month, $year, $actual_sunrise , $elevation) = ©fields; 


if ( ! defined ($elevation) ) { 

# only 4 fields were specified on the line - elevation not specified, so we use the last 

# specified elevation, exiting with an error if the latter does not exist 
if (! defined ($saved_elevation) ) { 

printf (STDERR "Input files do not specify elevation data on first data-containing line.\n"); 
exit 1; 

> 

$eIevation = $saved_elevation ; 

} else { 

# 5 fields specified on the line - save off the elevation 
$saved_elevation = $elevation; 


my $date = DateTime ->new( 
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day 

month 

year 


=> $day, 

=> $month, 
=> $year 


); 


# zero-based day offset from Winter Solstice 

my $day_offset = $date->delta_days($winter_solstice)->in_units(^daysO ; 

# various auigles 

my ($psi, $delta, $theta) ; 

# day of year angle 

$psi = ($day_offset / 365) * pi2; 

# calculate Pole Star elevation in radiains 
$delta = deg2rad($elevation) ; 

# sunrise angle 

$theta = deg2rad($actual_sunrise) ; 

my $sin_axial_tilt = -sin($theta) * cos($delta) / cos($psi); 
my $axial_tilt = rad2deg( asin($sin_axial_tilt) ); 

# exclude any estimates < 0 as abnormally inaccurate due to cos($psi) and sin($theta) becoming small, 

# thus introducing potentially large percentage errors, even though $psi and $theta are good 

# direction estimates. Inaccurate estimates on the +ve side will still be included though, 
if ($axial_tilt >= 0) { 

$total += $axial_tilt; 

$data_counter++; 

} 

printf ("7,6.2f7,5s7,5.2f7,9s7,5. If7,14s7,7.2f \n" , rad2deg($psi) , cos($psi), $actual_sunrise , $axial_tilt) 

> 

printf("7«d data points\n", $data_counter) ; 

printf ("Average axial tilt = 7»6.2f\n", $total/$data_counter ) ; 


Sample output : 


C:\>perl calc-axial-tilt.pl iedinburgh.txt 


PSI 

cos(PSI) 

Actual Sunrise 

Axial Tilt 

147.95 

-0.85 

39.0 

24.57 

157.81 

-0.93 

43.0 

24.36 

167.67 

-0.98 

46.0 

24.35 


Edinburgh 

37 data points 

Average axial tilt = 23.88 
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C.3 Calculate Analemma 


# Outputs a table of solar altitude and azimuth values for a given location on Earth and time of day. 

# The output can be used directly as input to TIKZ/pgfplots package to produce an analemma graph of altitude versus 

# azimuth (cuialemma-graph.tex). Uses a circular orbit/constaint speed model to calculate approximate orbital angles 

# for days of the year. 

# Usage : 

# perl calc-analeinma.pl <Latitude> <ClockTimeExDST> <LongitudeDiff> <DayStep> <NumberOfDays> ... 

# ... <AdjustAzimuthAngleRcLnge> <PerihelionDffset> 

# <Latitude> = latitude in degrees. 

# <ClockTimeExDST> = clock time in your time zone in decimal hours, excluding any DST shift. 24 hour clock. 

# <LongitudeDiff> = difference in degrees between your longitude euid longitude of the defining meridiem for your time zone. 

# The defining meridian for the time zone UTC+<n> is at longitude a multiple of <n> times 15 degrees (<n> +ve or -ve). 

# Your time zone^s clock (ie civil time, ex-DST) gives the mean solar time (MST) at your time zone’s defining 

# meridian longitude - the MST at your longitude at <ClockTimeExDST> is the latter plus (4 mins * <LongitudeDiff>), 

# where <LongitudeDiff> is +ve for easterly, amd -ve for westerly. Default value of <LongitudeDiff> is 0, in which 

# case the amalemma is calculated for the time zone’s defining meridian longitude rather tham your own longitude. 

# Ref : https://web.archive.Org/web/20190820154547/http://aa.usno.navy.mil/faq/docs/eqtime.php 

# <DayStep> = step between days (default value = 1). 

# <NumberOfDays> = length of period covered, starting from winter solstice, 365 for whole year (default value = 365). 

# <AdjustAzimuthAngleRange> = if set to true (1) adjust azimuth angle range from [0, 360) to (-180, 180] (default = 0), 

# use if angles congegrate around zero, to avoid discontinuity in the graph. 

# <PerihelionDffset> = offset in days of perihelion from winter solstice (default 14 = 2020 value). Used for sine wave 

# approximation of the eccentricity component of the Equation of Time. 

# Number of command line arguments = 2-7, the first two always required, and the remaining five optional in the 

# order specified. If ainy optional argument is specified as then it is set to the default value. 

use strict; 

use warnings ; 

use 5.012_003; 

use Math::Trig; 

use Math::Trig qw(:pi); 

my $AXIAL_TILT = deg2rad(23.44) ; 

if (OARGV == 0 $ARGV[0] eq "-h") { 
print ( "Usage:\n" . 

"perl calc-analemma.pl <Latitude> <ClockTimeExDST> <LongitudeDiff> <DayStep> <NumberDfDays> " . 
"<AdjustAzimuthAngleRange> <Perihelion0ffset>\n" . 

"First 2 parameters required, parameters 3-7 optional with default values : 0, 1, 365, 0 (false), 14\n" . 
"Specifying an optional parameter as ’.’ selects default value.\n" 

); 

exit ; 

> 

# number of command line arguments must be 2-7 
if (OARGV < 2 OARGV > 7) { 

print(STDERR "Invalid command line : 2-7 parameters required\n") ; 
exit 1; 

} 

# defaults 

my ($longitude_diff , $day_step, $num_days, $adjust_azimuth_angle_range , $perihelion_offset ) = (0, 1, 365, 0, 14); 

if (OARGV == 7) { 

$perihelion_offset = pop OARGV; 

$adjust_azimuth_angle_range = pop QARGV; 

$num_days = pop OARGV; 

$day_step = pop QARGV; 

$longitude_diff = pop @ARGV; 

} elsif (OARGV == 6) { 

$adjust_azimuth_angle_range = pop OARGV; 

$num_days = pop OARGV; 

$day_step = pop QARGV; 

$longitude_diff = pop OARGV ; 

} elsif (OARGV == 5) { 

$num_days = pop OARGV; 

$day_step = pop OARGV; 

$longitude_diff = pop @ARGV; 

} elsif (OARGV == 4) { 

$day_step = pop QARGV; 

$longitude_diff = pop OARGV ; 

} elsif (OARGV == 3) { 

$longitude_diff = pop QARGV; 

} 

# if ciny optional argument was specified as ’. ’ then reset it to default value 
if ($perihelion_offset eq ’.’) { 
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$perihelion_offset = 14; 

> 

if ($adjust_azimuth_aiigle_raiige eq { 

$adjust_azimuth_aiigle_raiige = 0; 

> 

if ($nuin_days eq ’.’) { 

$nuin_days = 365; 

> 

if ($day_step eq { 

$day_step = 1; 

} 

if ($longitude_diff eq { 

$longitude_diff = 0; 

} 

my $clock_time_ex_dst = pop OARGV; 

my $latitude_degrees = pop QARGV; 

# Calculate mean solar time corresponding to clock time, in decimal hours. 

my $mean_solar_time = $clock_time_ex_dst + ($longitude_diff / 15); 

# latitude in radiauis 

my $delta = deg2rad($latitude_degrees ); 

my ($sin_delta, $cos_delta) = (sin($delta) , cos ($delta) ); 


# sub sin_cos_lambda 

# ================== 

# Returns the sin/cos of the solar declination for a day of the year. 

# Usage : sin_cos_lambda($day_offset) 

# $day_offset = day offset from winter solstice (starting from 0) 
sub sin_cos_lambda { 

my $day_offset = $_ [0]; 

# day of year angle 

my $psi = ($day_offset / 365) * pi2; 

my $sin_lambda = -sin($AXIAL_TILT) * cos($psi); 

my $cos_lambda = sqrt(l - $sin_lambda**2) ; # note lambda is always in rainge [-90, 90] 
return ($sin_lambda , $cos_lambda) ; 

} 

# sub angle_from_sin_cos 

# ====================== 

# Returns an angle in the rainge [0, 360) from its sine and cosine values. 

# Usage : angle_from_sin_cos($sin, $cos) 
sub angle_from_sin_cos { 

my $sin = $_ [0]; 
my $cos_negative = ($_[!] < 0); 
my ($angle, $ql_angle, $mod_sin) ; 

$mod_sin = abs($sin); 

$ql_angle = asin($mod_sin) ; 

if ($sin >= 0) { 

if ($cos_negative) { 

$cLngle = pi - $ql_angle; 

} else { 

$aLngle = $ql_angle; 

} 

} else { 

if ($cos_negative) { 

$aLngle = pi + $ql_angle; 

} else { 

$cingle = pi2 - $ql_angle; 

} 

} 

return rad2deg($angle) ; 

} 

# sub solar_hour_angle 

# ==================== 

# Return the solar hour angle in radiauis using equation of time. 

# Usage : solar_hour_auigle($day_offset) 
sub solar_hour_angle { 

my $day_offset = $_ [0]; 
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# EOT component amplitudes in decimal hours 
my ($EC_amplitude , $OB_amplitude) = (0.1277, 0.1645); 


my $EC_component = $EC_amplitude * -sin( ($day_offset - $perihelion_offset ) * pi2 / 365 ) 
my $QB_component = $OB_amplitude * -sin( $day_offset * 2 * pi2 / 365 ); 

my $EOT_hours_adjustment = $EC_component + $OB_component ; 

# uncomment one of these to investigate effect on cinalemma of the separate EOT components 

# $EOT_hours_adjustment = $EC_component; 

# $EOT_hours_adjustment = $OB_component; 

# $EOT_hours_adjustment = 0; 

my $apparent_solar_time = $mean_solar_time + $EOT_hours_adjustment ; 
return ( ($apparent_solar_time - 12) / 24) * pi2; 

} 


# the hash allows the output to be used directly as a table in TIKZ/pgfplots 
printC'Day Altitude Azimuth\n"); 

print ("#- - - \n"); 

for(my $day_offset = 0; $day_offset < $num_days; $day_offset += $day_step) { 

# altitude 

my ($mu, $sin_mu, $cos_mu) ; 

# azimuth 

my ($phi, $sin_phi, $cos_phi) ; 

# solar declination 

my ($sin_lambda , $cos_lambda) = &sin_cos_lambda($day_offset) ; 

# solar hour aingle for clock time $mecLn_solar_time 
my $tau = &solar_hour_angle($day_offset) ; 

my ($sin_tau, $cos_tau) = (sin($tau), cos($tau)); 

$sin_mu = $cos_tau * $cos_delta * $cos_lambda + $sin_delta * $sin_lambda; 

$cos_mu = sqrtd - $sin_mu**2) ; # note mu is always in range [-90, 90] 

$mu = rad2deg( asin($sin_mu) ); 

# $cos_mu == 0 means sun is at zenith or nadir, therefore azimuth is undefined 
if ($cos_mu != 0) { 

$sin_phi = -$sin_tau * $cos_lambda / $cos_mu; 

$cos_phi = ($cos_delta * $sin_lambda - $cos_tau * $sin_delta * $cos_lambda) / $cos_mu 

} 

$phi = &angle_from_sin_cos($sin_phi , $cos_phi) ; 

if ($adjust_azimuth_angle_range && $phi > 180) { 

$phi -= 360; 

} 

if ($cos_mu != 0) { 

printf ("y,3d7,6s7,6.2fy,8sy,6.2f\n" , $day_offset, $mu, $phi); 

} else { 

printf ("y,3dy,6sy,6.2fy,8sy,s\n" , $day_offset, $mu, "Undefined"); 

> 

} 


Sample output : 

C:\>perl calc-analemma.pl 37.98 16 -6.27 I perl -S tee.pl analemma.dat 
Day Altitude Azimuth 


#— 


0 

1 


10.27 

10.35 


229.12 
229.04 


364 


10.21 


229.21 
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Analemma graph generator amalenuna-graph.tex : 


Xdocumentclass [a4paper]{article} 

\usepackage [margiii=0 . Sin] {geometry} 

\usepackage {tikz} 

\usepackage {pgf plots } 

\pgfplotsset{width=12cm,compat=l .14} 

\newcoinincLnd{\dg}{'‘{\circ}} 

\begin{docuinent} 

\begin{tikzpicture} [ 
node font=\Large 

] 

\definecolor{Grid Color }{HTML}{aaaaaa}; 

\definecolor{Analemma Color}{HTML}{f f 0000} 

\begin{axis} [ 
clip=f alse , 
axis equal =true , 
xlabel={Azimuth ($\dg$)}, 
xlabel style={anchor=north, below=8pt}, 
ylabel={Altitude ($\dg$)}, 
ylabel style={anchor=south, above=8pt}, 
grid=both, 

grid style={Grid Color, line width=0.3pt, opacity=0.3}, 

] 

\addplot [ 

Analemma Color, 
only marks , 
mark size=l.lpt, 

] tableCx index=2, y index=l, skip first n=0] {cinalemma.dat} — cycle 
\path (axis description cs:0.5, 1.1) node [] {$<$Title$>$}; 

\end{axis} 

\end{tikzpicture} 

\end{document } 
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Appendix D 


Rodrigues Rotation Formula 

In the following EV denotes the set of all Euclidean vectors. Rotations about an axis follow the right hand rule convention. 


Theorem. Given a unit vector u G EV then the rotation r' of vector r G EV about the axis u is given by : 

r' = (cos 0) r + (sin 0) (u x r) + (1 — cos 9) (u • r) u 
Proof. As shown in Figure D.l express r as a sum of components parallel and perpendicular to u : 


(37) 


r r r 

where r 11 = (u • l) u , 

and r = r —r,, = r—(u-r)u. 

Since we are using the right hand rule convention the image of rounder the rotation is given by : 


r^ = (cos 0) r^ + (sin 0) (u x r^) 

= (cos 9)y_ + (sin 0) (u x r) , since u x 


and clearly r' = r thus : 


= (u • r) u + (cos 0) (r — (u • r) u) + (sin 0) (u x r) 


which readily rearranges to (37). 


QED 





Figure D.l; Rotating Vector r by d About Axis u 
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Appendix E 


Latitude as Elevation of Pole Star 

In Figure E.l, an observer at latitude 6 above the equator sees an angle of elevation 7 of the Pole Star above their horizon. 
Because the distance OP is vastly much greater than the planet radius R the lines PA and PN are close to parallel, 
thus the angles of incidence 6 and 7 of the horizon line onto these two lines are close to equal. (And as we approach the 
equator the Pole Star falls to the horizon with 7 < i5, i5 —>■ sin “^(^)+,7 —>• 0+, and R <C OP). Thus 7 always well 
approximates d as a direction. 



Figure E.l: Latitude as Elevation of Pole Star 
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